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ABSTRACT 



To constrain the Galactic gravitational potential near the Sun (~1.5 kpc), 
we derive and model the spatial and velocity distribution for a sample of 9000 
K-dwarfs that have spectra from SDSS/SEGUE, which yield radial velocities and 
abundances ([Fe/H] & [a/Fe]). We first derive the spatial density distribution 
for stars of three abundance-selected sub-populations by accounting for the sur- 
vey's selection function. The vertical profile of these sub-populations are simple 
exponentials and their vertical dispersion profile is nearly isothermal. To model 
these data, we apply the 'vertical' Jeans Equation, which relates the observable 
tracer number density and vertical velocity dispersion to the gravitational poten- 
tial or vertical force. We explore a number of functional forms for the vertical 
force law, and fit the dispersion and density profiles of all abundance selected 
sub-populations simultaneously in the same potential, and explore all parameter 
co-variances using MCMC. Our fits constrain a disk mass scale height < 300 pc 
and the total surface mass density to be 67 ±6 M Q pc -2 at \z\ = 1.0 kpc of which 
the contribution from all stars is 42 ±5 M Q pc -2 (presuming a contribution from 
cold gas of 13 M pc~ 2 ). We find significant constraints on the local dark mat- 
ter density of 0.0065 ± 0.0023 M pc~ 3 (0.25 ± 0.09 GeVcm~ 3 ). Together with 
recent experiments this firms up the best estimate of 0.0075 ± 0.0021 M pc~ 3 
(0.28 ±0.08 GeVcm" 3 ), consistent with global fits of approximately round dark 
matter halos to kinematic data in the outskirts of the Galaxy. 

Subject headings: 
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Introduction 



Since the works of lOortl (Il932l . Il960l ). determining the Galactic gravitational potential 
from the distribution and kinematics of stars has been one of the most instructive problem 
in Galactic disk study. Comparing it with the mass density distribution of visible material 
(stars and gas), one can derive the distribution of dark matter in the Galactic disk. 



Comp arisons of such local dark matter estimations with rotation curve (IWeber fc de Boer 
20101 ) constrain the shape of the dark matter distribution. 

Q 



Kuiiken fc Gilmord f ll989aU b 



hereafter KG89a,b,c) pioneered a practical approach 



to constrain the Galactic potential near the Sun. They used K-dwarfs to estimate the 
surface mass density of the total gravitating mass and the mass density of dark matter by 
solving the 'vertical' Jeans' Equation (i.e. in the ID equation in the z-direction) , which 
correlates the space density and velocity of tracer stars. They found the identified (i.e. 
stellar and gaseous) surface density and the total surface mass density of all gravitating 
matter within \z\ < 1.1 kpc from the Gal actic plane near the Sun t o be 48 ± 8 M Q pc~ 2 and 



71 ± 6 M pc" 2 , respectively (KG89a,b,c 



Kuijken fc Gilmore 



1991 



hereafter KG91). The 



uncertainties are mainly caused by the measured errors in distance and velocities of the 
K-dwarfs and the technique of recovering the vertical force K z . 

After KG91 there have been various other determinations of the disc surface mass 
density based on different stellar tracer: iFlynn &: Fuchsl (1l994j . hereafter FF94) derived 



a sur face mass density of 52 ± 13 M & pc 2 for the known disc matter. With Hipparcos 



data, 



Korchagin et al. 



Holmberg fc Flynnl (120041 . hereafter HF04) and 



jgOOa l focused S 



,z|<350 pc 



42 ± 6 Mg, pc 



Bienayme et al. 



, by using red giants; 



(120061 ) used K giants 



and red clump stars to estimate the disk surface density again, and the values are 



J |z|<l.l kpc 



74 ± 6 M pc~ 2 and Xm<i.i 



kpc 



57 — 79 ± 6 M Q pc 2 , respectively. The 



rather broad range in values implies that better data and techniques are needed to improve 
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the measurement of the density distribution of Galactic disk. 

In this paper, using the large sample of K-dwarfs observed in SDSS/SEGUE 



( lYanny et al. 



20091 ). we re-determine the disc surface mass density and mass density of dark 
matter. In a number of aspects, this work follows KG89b, in particular in constraining the 
"vertical force" K z . However, the paper presents a number of new elements over previous 
studies in addition to the new large data set. We split the data into abundance selected 
sub-samples, which provide distinct probes of the same potential; and we simultaneously fit 
densities and kinematics using Markov-Chain Monte Carlo (hereafter MCMC) approach. 
Besides, we explore how the results depend on the functional forms for K z . 

This paper is organized as follows. In § [2, we describe the SEGUE/SDSS data, in 
particular we emphasize the spatial selection function, the distance determination, and the 
sub-samples of similar [Fe/H] and [cn/Fe]. In §[3] we lay out how we determine the tracer 
number density and vertical velocity dispersion, and practicalities of solving the Jeans 
equation. The fitting results are shown and discussed in § H] and § [5j with a summary in § [6j 



DATA 



For our analysis we aim for a sample of kinematics tracers that can be found in 



disk populations of various ages and metallicities, and which cove r distance o 



to 2 kpc from the Sun. In the Sloan Digital Sky Survey (SDSS; 



York et al 



0.2 kpc 



2000), 



K-dwarfs best satisfy these criteria. The Sloan Extens ion for Galactic Under standing 



20091 ) which 



and Exploration (SEGUE) is a subsurvey of SDSS-II (lAbazajian et al. 
operated from 2005 August to 2008 July, to probe the formation and evolution of our 
Galaxy. I t obtained uqriz imaging of some 3500 deg 2 of sky outside of the SDSS-I 



footprint (IFukugita et al. 



1996 



Gunn et al. 



1998 



2006 



York et al 



200ol : 



Smith et al 
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2002; 


Stoug 


iton t 


it al. 


2002; 


Abazaiian et al 


. 2003. 


2004. 


200 


5 


2009; 


Pier et al. 


2003 




Ivezic et al. 


2004; 


Adelman-McCarthv et al. 


2006 


> 


2007 




2008 




Tucker et al. 


2006 


), with 



special attention being given to scans of lower Galactic latitudes (|6| < 35°) in order to 
better probe the disk/halo interface of the Milky Way. Overall, SEGUE obtained some 
240,000 medium- resolution (R ~ 2000) spectra of stars in the Galaxy, select ed to explore 



the nature of stellar populations from 0.3 kpc to 100 kpc (lYanny et al. 



20091). The seventh 



data release (DR 7) is the final public data release from SDSS-II occurring in October 2008. 
SDSS-III, which is presently underway, has already completed the sub-survey SEGUE-II, 
an extension intended to obtain an additional sample of over 120,000 spectra for distant 
stars that are likely to be members of the outer-halo population of the Galaxy. Data 
from SEGUE-II hay e been distributed as part of the eighth public data release (DR 8) 



(Aihara et al. 



2011al ). The SEGUE Stellar Parameter Pipeline processes the wavelength- 



and flux- calibrated spectra gen erated by the standard SDSS spectroscopic reduction 



pipeline (jStoughton et al 



20021 ) . obtains equivalent widths and/or line indices for more 



than 80 atomic or molecular absorption lines, a nd estimates T P ff, lo g q, and [Fe/H] throu gh 



the application of a numbe r of approaches (see 



Lee et al. 



2008a 



2008 



Smolinski et al. 



201l|). 



b; 



Allende Prieto et al. 



2.1. K-dwarf selection 



KG89b suggested that the tracer star in the present context should have the following 
properties: 1) they are phase-mixed and in dynamical equilibrium in Galactic potential; 
2) they are common and distributed through all Galactic disc components; 3) the distance 
of tracer star can be well determined. Therefore, considering present SDSS/SEGUE 
observations, K-dwarfs are ideal stars which are used as tracer to measure the total disc 
mass. Compare to the SEGUE G-dwarfs (e.g. 



Bovy et al. 



2012al ). they have the advantage 
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that their minimal distance to be included in SDSS/SEGUE is almost two times closer. 



Yanny et al. 



(120091 ) list the following color and magnitude cuts for identification and 



targeting of K-dwarf candidates: 

14.5 < m r < 19.0, 
0.55 < (g -r) < 0.75. 

In order to get reliable estimates for the 'vertical' (z) component of each star's velocity, 
we added two more criteria, the existence of a good proper motion measurement and the 
existence of spectra in the data base 

error of proper motion > 0., 
S/N > 15. 



These criteria return 10,925 candidates in DR 8. However, DR 8 does not list [a/Fe], which 
is very important in subsequent analysis for K-dwarf candidates. Besides, proper motions 



of DR8 above declination 5 > 4 



calibration in DR8 (Aihara et al. 



° are worse than those of DR7 because of the astrometric 



2011bl ). Thus, we use the [a/Fe] and proper motions from 
DR 7 for all of these candidates. To eliminate K-giants, only stars whose logg > 4.0 are 
taken in the present work, which in the end leaves 9,157 stars. Fig. [1] shows the number 
density distribution of sample stars in [a/Fe] vs. [Fe/H] space. 



2.2. Abundance-selected sub-samples 



As shown in Fig. (TJ the distribution of the K-dwarf sample in the [a/Fe] — [Fe/H] 
plane is bi-moda L i.e. a m e tal-ri ch, q -deficient populatio n and a metal-poor, a-enhanced 
one, as found by lLee et al.l (120111 ) and iBovy et al.l (I2012al ) for SDSS G-dwarfs. 



Liu van de Venl (120121 ) and iBovy et al.l (l2012bl ) have shown that the kinematics of 



'mono'-abundance sub-populations are relatively simple. Therefore, we will split our sample 
of K-dwarfs into subdsets that are abundance-selected in the [a/Fe] — [Fe/H] plane (see 
Fig. [TJ black boxes) 

• metal-rich: [Fe/H] G [-0.5, 0.3], [a/Fe] G [0., 0.15] 

• intermediate metallicity: [Fe/H] G [-1.0, -0.3], [a/Fe] G [0.15, 0.25] 

• metal-poor: [Fe/H] G [-1.5, -0.5], [a/Fe] G [0.25, 0.50] 

These sub-samples contain 3672, 1416, and 2001 stars, respectively. Of course, all these 
sub-population stars move in the same gravitational potential. We will present the analysis, 
separating the three sub-populations. 



3. The Galactic potential near the Sun 



For the present context, we only consider the simple problem of solving the Jeans 
equation for a one-dimensional slab and then to determine gravitational potential and hence 
density of matter near the Sun (\R — R & \ = 1.0 kpc and \z\ < 1.5 kpc). Following KG89b, 
we assume that: the Galactic potential is the sum of contributions by observable matter 
and by dark matter and we assume the effective dark matter density is a constant. Hence, 



the "vertical" Jeans equation can be written as (IBinney &: Tremaine 



2008|) 



_d 



2 ] 



■0; 



dz 



where R@) is the vertical gravitational potential in the Sun's vicinity, i.e., we solve this 
equation for R = constant. Here, v^z) is the vertical number density of the tracer stellar 
population, and u z {z) the vertical velocity dispersion of exactly those tracers. 



The first derivative of $(z,Rq) is the gravitational force perpendicular to the galactic 
plane K z (z) 

K.fr^-WZM (2) 
az 

The (ID) Poisson equation relates the potential to the disk vertical mass density, 

4:irGp tot {z,R Q ) = — (3) 

where ptot{ z , R&) is the total mass density of all luminous components and the dark matter 
contribution. 

Combining Eq. (TJ Eq. |2l and Eq. |3l one obtains 

A [^(z)a z (z) 2 ] = -v*(z)K z (z,R Q ) (4) 

A-nGp tot (z } R Q ) = ~K z {z, R & ) (5) 

To solve Eq. |4] for a given tracer sub-population and derive the local mass density 
distribution by using Eq. [51 we proceed as follows: 



i) We determine v*(z) and cr z (z) for each of the three tracer populations directly from 
the observations. 



ii) We presume the v±{z) is a simple exponential ( IBovy et al.ll2012af ) of unknown scale 



height h z . Here h z = [h Zj i,h Zj 2,h z ^] correspond to the metal-rich, inter-intermediate 
metallicity, and metal-poor sub-samples, respectively. 

iii) We adopt a parameterized form of K z (z, R & \ p). 

iv) For each [p, h z ], we fit for the best u*(z) and cr z (z), and get the probability of the data 
given the parameters, i.e., V(a z , v± \ p, h) 

v) We use a MCMC technique to sample the likelihood of the parameters and find the 
best [p, h z ] . 
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During the whole calculation, we assume a parameterized K z model instead of 
gravitational model, because it makes the calculation easier. After getting a reasonable K z 
form, gravitational potential can be derived by integration of K z and mass density can be 
determined through the first derivative of K z . Here we will describe the details step by step. 



3.1. Tracer number density 

3.1.1. Distance estimates and coordinate system 

To obtain z and v z of our sample stars in the Solar vicinity, we need to estimate 
their distance first. Given the color indices (g — r and g — i), the r-band magnitude, and 
metallicity, the absolute magnitude is estimated by fitting to fiducial color-magnitude 



relations, calibrated through star cluster spectroscopy. The fiducial s equences 



and (g — i, r), based on YREC+MARCS isochrones, as described in 



AnetaL 



or (q — r, r) 



(|2009|), are 



adopted. Then it is straightforward to determine the distance D = 10 s r+1 pc. 

Such distances were estimated by fitting the fiducials for the two colors above 
separately, yielding an average difference and standard deviation among these two distance 
estimates of 0.006 ± 0.077. We use the average distance and the standard deviation from 
color-magnitude (hereafter CM) diagrams of (g — r, r) and (g — i, r) as the mean distance 
and its error, respectively. 

With D estimated as above and observed Galactic longitude and latitude(Z, b), we get 
the stars' position in Galactic cylindrical coordinates. We adopt a Galactocentric cylindrical 
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coordinate system in which z-axis towards the North Galactic Pole (b ^ 90°), that is, 



R = 

9 = arctan 
Z = Dsm{b) 



D 2 cos 2 (b) - 2DR & cos(6) cos(Z) + R% 
D cos(6) cos(Z) — Rq 



D cos(fe) sin(Z) 



(6) 



19931 ). In the 



where Rq = 8.0 kpc is the adopted distance to the Galactic center ( jReidl 
following works, only the sample stars in the range 0.3 kpc < \z\ < 1.2 kpc are used. This 
z-range is limited by the bright limit of SDSS data at the low end and by the legitimacy of 
treating the problem as 1-D (z-direction) for large \z\. 



3.1.2. Selection function 

Any dynamical analysis need accurate knowledge of the spatial distribution of the 
kinematic tracer. Therefore, we need not only a well-defined spectroscopic sample of 
K-dwarfs with known fluxes and distances (or he ights above the m id-plane) estimates, we 



2012a ). which is what we lay 



also need to understand their selection function fjBovy et al. 
out there. 

In this coordinate system, each sub-sample can be divided in eight bins according to R 
and AR = R m+ i — R m = 500 pc. In each AR range, R is considered as a constant, and each 
sub-sample is also divided in 12 bins in z direction. The width of the z bin, Az = Zi+i — Zi, 
is Az = 100 pc. Here m and i are indices of R and z bins, respectively. In each z-R box 

u ^b (R \ N suh (Rm,Zi) 

where N suh (R m , Zi) and V^- h (R m , z{) are respectively the star number and the effective 
volume in a (zi,R m ) box of a certain sub-sample. Because every targeted star of a given 
(r, g — r) has a possibility to be in one particular sub-sample, each line of sight is a part of 
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the search volume for every abundance sub-sample. In our effective volume calculation, all 
lines-of-sights of all K-dwarfs are included, i.e. 

v:i h (R m , Zi ) = J2 V ^ ( 8 ) 

9=1 

where n su b is the total number of lines-of-sights of a given sub-sample and V e ff, q is the search 
volume of each line-of-sight. 

It is in the calculation of the effective volume that the selection function enters 
explicitly. Some aspects of the SEGUE selection function are obvious. The apparent 
magnitude range brackets the possible distances of K-dwarfs, and the most nearby stars 
will be preferentially the brightest and coldest stars. At a given m r and g — r, the effective 
survey volume is smaller for lower metallicity (hence less luminous) stars. Moreover, the 
SDSS targeting strategy implies that stars at lower latitude have a smaller probability of 
ending up in the spectroscopic sample, since a fixed number of targets is observed along 
each line of sight. In general, we need to derive the selection function, which gives the 
probability that a star of given M r , g — r, [Fe/H] ends up in the samples, as a function of 
D and (/, b) 

Each SEGUE pointing which corresponds to an angular on-the-sky radius of 1.49°is 
observed with two plates, a SEGUE bright plate that ta rgets stars with 14 .0 < m r < 17.8 



20091 ). In each plate, 



and faint plate that targets stars with 17.8 < m r < 20.1 ( lYanny et al.l 
the distribution of spectroscopic sample in a CM box that satisfies the selection criteria is 
^spectj, g — r), while the distribution of all photometric stars within the same plate that 
satisfy the same CM cuts is n p ^ oto (r,g — r). Integration of these distribution over all CM 
cuts results in A^ spec and N photo , leading to a plate weight 

Opiate, = (9) 
* ' photo j 

where j is the index of plate. This plate weight should be a part of selection function. 
Fig. [2] presents the cumulative distributions of spectroscopic and photometric samples in 
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the same selection criteria. It is clear that the selection weight for faint plates has an 
apparent magnitude dependence because of the dependence on signal-to-noise ratio of the 
probability of obtaining a good spectrum but no strong color dependence. Therefore, we 
make an approximation for this selection function, that is, it is a function of m r . 

This function is defined for all distance moduli, which is one for distance moludi that 
correspond to apparent magnitude insides of plates' magnitude range and zero for other 
situation. We combine the selection function in to the calculation of V e g, q , that is, the 
average metallicity of one particular sub-sample <[Fe/H]> is adopted to estimate the 
possible absolute magnitude, M q , distance, D q , and height above mid-plane, z q , of a given 
(r, g — r). Then 

V eS , q = | • W plate , q ■ W q • |^ • i[cat(&, -9)- cot(6 9 + 6)} ■ (z 3 nppcr - zf oweT ) (10) 



wher e 6 = 1.49°is the diameter of SEGUE (and SDSS) spectroscopic plate ( lYanny et al 



2009|). Besides, 

1, for Zi < Z q < Z i+ i 

0, otherwise 



W q = { ' (n) 



■^upper 

= Min(z i+ i, -D r , m ax * sin(6,)) (12) 
z lowev = Max(zi,D rMn *sm(b q )) (13) 



D r , min = lO^f^^pc (14) 

A,max = 10^^ + 1 P C (15) 

where r min and r max are the magnitude limits of each plate. For SEGUE bright plate, 
7 >in = 14.0 and rvn ay = 17.8, and for SEGUE faint plate, r min = 17.8 and r max = 20.1 



(jYanny et al 



20091 ). In each Z{ bin, the error bar of the tracer number density arises from 
the star co unt Poisson varia nce and it is estimated by means of Monte Carlo bootstrap ping 
(§ 15.6 of IPress et al.ll2007l ). 
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Fig. [3] presents the resulting (R,z) map of the tracer number density for each 
sub-population. From top to bottom, those plots correspond metal-rich, intermediate 
metallicity, and metal-poor sub-sample; from left to right, plots represent effective star 
number, effective volume, and natural logarithm of number density in each z-R box. In 
this figure, one can see that the scale-height of metal-rich sub-sample is shorter, but the 
scale-length of the sa me sub-population is longer. This result is similar as the studies of 



Galactic structure of 



Bovy et al 



( )2012d ). If this number density result is projected onto z 



direction, we get the vertical profile of the tracer number density. 



3.2. Vertical velocity dispersions 



To estimate cr z (z, [Fe/H], [a/H]), one could simply calculate the standard deviation 
of the observed v Zy i values for each z bin of each sub-population. However, this would 
encompases the contribution arising from t he measurement e r rors. Therefore, we use 
maximum likelihood technique described in 



van de Ven et al. 



(120061 . Appendix A) to 



estimate the intrinsic velocity dispersion, corrected for all individual velocity errors. 

In each z bin, the intrinsic velocity distribution of the stars is assumed as C(v z ). Each 
stellar velocity v Z) i in this bin is the product of C(v z ) convolved with a delta function which 
is broadened by the observed uncertainties 8v Z) i. For all N stars in the bin, the likelihood 
can be defined by 

»+oo 



N 



1 / Vx,i—"z 



C(v) 



2n<J; 



(16) 



where S£ is a function of mean velocity v z , velocity dispersion o z and possible higher-order 
velocity moments. The velocity distribution C(v z ) can be recovered by maximizing the 
likelihood £ £ '. To describe the velocity distribution, we u se a parameterized Gauss-Hermite 



(GH) series (Ivan der Marel fc Franx 



1993; 



Gerhard 



19931 ) because it has the advantage that 
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it only requires the storage of the velocity moments (v z , a z , hs iZ , h^ z , . . .) instead of the full 
velocity distribution and the convolution of Eq. [TB] can be carried out analytically. This 
makes it feasible to apply the method to a large number of discrete measurements and to 



estimate the uncertainties on the extracted ve 
bootstrap method (§ 15.6 of 



Press et al 



2003). 



ocity moments by means of the Monte Carlo 



Fig. H] shows an example of a vertical velocity dispersion fit in 600 pc < \z\ < 700 pc for 
the three sub-populations. This figure shows comparisons of different fitting methods. Red 
solid curves are the results of Gaussian fits without considering observed errors. Green and 
blue curves represent three and four moments GH fits using the method described above. 
For all \z\ bins of all sub-populations, the vertical velocity dispersion is overestimated by 2 
- 4 km/s without considering observed errors. 

For the metal-rich, intermediate metallicity, and metal-poor sub-samples, the measured 
o"2,obs(z) are in the ranges 17 - 26 km/s, 25 - 37 km/s, and 35 - 42 km/s, showing a 
increase with \z\. The u z ,obs{z) of different sub-populations are slightly increasing with the 
increasing of \z\. Moreover, 0" ZiO b s (-2) also increases with the decline of metallicity and the 



increasing of |a /Fe| (see Fig. El fj f 



ed ci rcles). We compare our observed cr zohs (z) with 



recent study of iLiu fc van de Venl (120121 ). in which they adopted SEGUE G-dwarf of the 



Solar neighborhood. They used the splits at [a/Fe]=0.25 and [Fe/H]=-0.6 to divide their 
sample into three groups: group I of a-young metal-rich stars which correspond to our 
metal-rich and intermediate metallicity sub-samples, group II of a-old metal-rich stars, and 
group HI of a-old meta l-poor stars which correspond to our metal-poor sub-sample (Fig. 3, 



Liu fc van de Yen 



20121 ). For the three groups, they found that 



;he oy nhs(-z) is about 20 



40 km/s, 35 - 50 km/s, and 40 - 55 km/s, respectively (Fig. 4, 



Liu fc van de Ven 



These cr Zfi b s (z) are in good agreement with our measurements. iBovy et al. 



(12012c 



20121 ). 



) also 



studied the vertical motions of mono-abundance sub-populations in the Galactic disk by 
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G-dwarf, and found that 0" 2iO bs(^) varies from 15 — 25 km/s, 20 — 30 km/s, and 30 — 
50 km/s for stars in our metal-rich, inter-intermediate metallicity, and metal-poor selection 
bins, respectively, which are consistent with our measurements. 



3.3. Functional forms for K z 

Given the tracer number density and vertical velocity dispersion probes, we investigate 
different parameterized K z (z) forms to solve Eq. [T]and derive the corresponding expression 
for cr z (z). K z (z) is contributed both by the baryonic components and the dark matter. 
According to the discussion of KG89b, all K z (z) forms should satisfy the following 
assumption: (1) the halo and disk contributions are degenerate when \z\ <C scale height; 
(2) disk contribution is approximate from a thin mass sheet when \z\ ^> scale height. We 
considered five different parameterized forms of K z (z), two of which we discuss below, while 
the other three are given in the Appendix, 



i) KG89 model 



y z + y 

* r^~, — 2 s as 

v z + 4 



+ 47rGp DM z (17) 



K z = 2nG 
ii) Exponential model 

1 - exp(-— ) + E gas I + ATrGp DM z ( 
z h J J 

Heres E* is the surface mass density in stars, E gas = 13.2 M pc~ 2 is the s urface mass 



K 7 = 2txG { E 



densit y of gaseous components which is taken from the disc mass model of 



Flynn et al. 



( 20061 ). and pdm is the effective dark matter contribution, and Zh is the scale height of mass 



density. 

In these models we use E* + E gas rather than pdisk, because pdisk and Zh are degenerate, 
i.e., one can not constrain p<jisk and Zh separately at the same time. 
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For the exponential model, the corresponding vertical dispersion profile can be derived 
analytically: 



1 z h . z , 
1 — - — exp 



4nGp DM h™ h (z + h 



,sub\ 



For the other K z forms we infer <7 2)SU b through numerical integration. 



3.4. Obtaining the PDFs for the model parameters 



(19) 



The predictions for <J zsuh (z) are based on 6 parameters for each K z (z, R Q ) model, which 
may have considerable covariance. At the same time, we have a parameterized model for the 
vertical tracer density profile, and we want to derive constraints on the gravitational force, 
marginalized over the scale height parameters. To explore the parameter space efficiently 
and to account for parameter degeneracies, we use a MCMC approach, which provides a 
straightforward and powerful way of estimating the probability distributions (PDF ) for all 
the model parameters and their degree of uncertainty (§ 15.8 of 



Press et al. 



20071 ). 



The objective is to calculate the probability of the data given the parameters, where the 
'data' are the tracer number density and the vertical velocity dispersion profiles described 
above. The probabilities for u z {z) and v*{z) in each \z\ bin are independent, therefore, the 
logarithm of the likelihood of the data given the parameters can be written as 



N 



N 1 

El 



- 1 



fat,- 



^obs, 



-Voba,3 



(20) 



where N is the total number of \z\ bin, p are the parameters of gravitational model, and 
h = [h Zt i, h Zi 2, h z3 ] are the scale heights of tracer number density of different sub-populations. 



-17- 



The trace number density profile is modeled by single exponential, 



v*, a (z) = Vq, s x exp y~—J ( 21 ) 

where the subscript s denotes the different abundance-dependent sub-populations. In order 

to reduce the number of parameters by eliminating z/ 0jS from Eq. [2U we calculate uq s for 
each sub-population by 

-jr- 1 = 0, (22) 
ov , a 

yielding 



N , v / JV 



^ = E^,/E^^(f] ■ (23) 

All other parameters are estimated by running a MCMC chain of typically 50,000 
steps, after a burn-in. In the Markov chain, the parameters of the gravitational model and 
scale heights of tracer densities [p n , h z>n ] are chosen from a multivariate normal distribution 
of [p n -i, h n _i], with the covariance of [p, h z ] chosen to yield a total acceptance rate in the 
MCMC Chain of about 20% ~ 30%. We use the results of all steps to represent [p, h z ] and 
their distributions, by counting the rejected steps as repeat ones. After sampling the PDF 
for parameters for each K z model, Xtot between the observed values (cr 2j0 b s and z>v |0 bs) an d 
model predicted ones (a 2imo d and z4, mo d) is calculated to explore how the models form the 
different K z families compare: 



2 2 | 2 



s=l 



N / \ 2 N , x 2" 

>y ^ / &z,j,s,mod &z,j,s,obs j ^ ^ > / ^j'^mod ^j,s,obs 

j V Z ,j ,S ,obs J j ^^J; s ;°bs 



(24) 



where s is the indices of sub-sample, j is the number of \z\ bins. 
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4. Results 

The above procedure results in PDFs for the K z and v+ parameters, which are presented 
in Figs. [5] — [TTJJ The best fitting results and recovered parameters for all K z gravitational 
force law models are listed in Table [TJ 



4.1. Vertical tracer density profiles 



Fig. [5] shows the fits to the vertical tracer density profiles. Filled circles are 
values derived directly from observations, which are the average value of R columns 
(7 kpc ^ R ^ 9 kpc) of the right panel of Fig. [3j Dashed lines are the predictions for u+ )S (z) 
(Eq. l2Tj) for two different K z models, and the shaded regions represent the samplings of 
the last 200 steps in the MCMC chain. Red, green, and blue symbols represent metal-rich, 
intermediate metallicity, and metal-poor sub-population, respectively. Each sub-population 
has a simple and single exponen tial structure, and the scale height is incre asing with the 



decline of metallicity (as seen in 



Liu van de Yen 



2012 



Bovy et al 



2012cf ) 



The joint parameter PDFs show that the tracer scale height, h z is independent of the 
K z model parameters (see Fig. [9]). The scale-heights we find for different abundance-select 
tracer sub-populations are 230 - 260 pc, 450 - 510 pc, and 800 - 1000 pc (see Table [TJ, for 
all K z families. The minor differences caused by different models are within the 68% central 
region of the h z PDF. 



These finding s are in accord with recent results by 



Bovy et al 



(]2012ch and 



Liu fc van de Venl ( 120121 ). In their work, they adopted G-dwarfs from SDSS/DR7 to 
analysis the relationship between the distribution of Galactic disc stars and their metallicity 
in each narrow [a/Fe]-[Fe/H] box. They concluded that the scale-height increases 



continually from 200 pc to 1200 pc with the decline of metallicity. 
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4.2. Vertical dispersion profiles 



As discussed in I3.2[ the a z , bs(z) shown in Figure 6, have already been corrected for the 
velocity errors of the individual stars within each \z\ bin. In these plots, we also summarize 
&z,mod( z ) predicted by different K z models. The results reflect simultaneous fit to all three 
sub-populations. If the (T z ,obs( z ) is fitted for each sub-population separately, there is some 
slight tension between the fits to the individual sub-populations and the simultaneous fit. 
The difference of <T 2 ,mod(^) given by different models is smaller than the error of observed 
Cz,dbs{z), as in the case of the scale-height predictions of vertical number density. That is, 
for the three sub-population, a Zjmo d(z) is in the range of 16 - 20 km/s, 25 - 29 km/s, and 35 
- 43 km/s (see Table [TJ, respectively. 



W e note again, that pr evious studies of K 7 [z] by KG89a,b,c, HF04 



(|2003|) 



Garbari et al.l (120121 . hereafter G12) and 



Bovy fc Tremaind (120121 . hereafter BT12) 



2. 



Siebert et al. 



had not, or could not, split their tracer samples a priori into abundance sub-bins with nearly 
isothermal o~ z (z) profiles. In samples with a wide metallicity range o~ z (z) will rise with 
\z\ simply because the mix of mono-abundance populations will change. To disentangle a 
rising u z (|z|) due to a dark matter halo from the effect of population mixing, these studies 
had to model the metallicity distribution, which turned out to be an important source of 
systematic errors. Fitting simultaneously to populations that have been split a priori by 
their abundances, as done here, reduces this source of systematics. 



4.3. The vertical gravitational force K z 

Fig. [7] shows K z (z) implied by our fits using the KG89b and exponential families for 
the 'vertical force law'. The dashed fat lines show K z (z) for the most likely parameters, 
and the band of grey lines show a la sampling of the PDF for K z (z). For both the KG89b 



-20 - 



and the exponential families of K z (z) show a number of generic features in Fig. [7J they 
start out with a finite value for small \z\, as we have fixed a prior contribution from a thin 
layer of the cold gas with ~ 13 M pc~ 2 ; then K z rises steeply to ~ 300 pc (reflecting the 
mass scale height of the stellar disk) and then flattens out; beyond ~ 500 pc the slower rise 
in K z reflects the dark matter halo term. 

One the basis of K z (z), we can derive the total surface density at \z\ = 1.0 kpc. We get 
£S 9 |<i.o kpc = 67 ± 6 M pc" 2 for the KG89 model and Egj,,^ kpc = 66 ± 8 M p C - 2 
for the exponential model. It is clear that Stot,b|<i.o k pc is robustly constrained, irrespective 
of the K z model. Our modelling also constrains the mass scale height, 245^245 P c f° r 
KG89 and 20013,0°, f° r the exponen tial model ( see F ig. ED. Within these uncertainties, this 



is consistent with the 180 pc from 



Hill et al. 



390+330 pc of 



Siebert et al. 



(J1979J) by analyzing A and F dwarfs and 



(120031 ) . found by using high resolution spectral data of red 
clump stars. Taken together this forms up the picture that there is a dominant mass layer 
near the disk mid-plane (presumably baryonic) that is rather flat. 

For comparison with literature studies, we also explicitly estimate the total surface 
density of baryonic matter within \z\ < 1.1 kpc, £bar,|«|<i.i kpc — + S gas , which we find 
to be 5™ <L1 kpc = 55 ± 5 M p C - 2 , and E^ p |z|<11 kpc = 54 ± 8 M p C - 2 , respectively. 
These values are slightly higher than 48 ± 8 M pc~ 2 derived by KG89b, but are in perfect 
agreement with the value of 52 ± 13 M pc~ 2 of FF94 and 53 ± 6 M pc~ 2 of HF04. 

The local mass densities of dark matter recovered by our models are Pdm 89 = 
0.0065 ± 0.0023 M pc" 3 (0.25 ± 0.09 GeVcm~ 3 £ and p^M = °- 0060 ± °- 0020 ^©pc~ 3 
(0.23 ± 0.08 GeVcm -3 ), marginalized over the other parameters.. We explored the 
parameter degeneracies by MCMC (see Fig. [TOl) . which recover the known degeneracy 



x l GeVcm" 3 ~ 0.0263158 M pc" 
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between surface mass density of stars and volume mass density of dark matter: lower E* 
correlates with higher Pdm, at an approximately constant total surface mass density 



5. Discussion 



We have presented an analysis of the vertical Galactic potential at the Solar radius, 
drawing on ~ 9000 K-dwarfs from SDSS/SEGUE. In many ways, the analysis followed 
the approac h initially laid out by KG89 and implemented by several other groups in the 



mean-time (ISiebert et al. 



2003 



HF04, G12, BT12). In comparison to most previous studies, 
our analysis has a number of new elements: we have a substantially larger sample than 
previous analyses; we have taken explicit account of the abundance-dependent spatial 
selection function of our sample for the analysis; we have simultaneously fit for several 
abundance-selected nearly-isothermal sub-populations that 'feel' the same gravitational 
potential, and we have matched the kinematics and the spatial distribution simultaneously; 
we have explored to which extent the choice of the functional form for K z affects the results. 

As laid out in the previous Section, we find good constraints on E to t,|z|<i.o k pc an d 
some constrains on thickness of the disk mass layer: Zh- £tot,|«|<i.o k pc = 67 ± 6 M pc~ 2 
and Zh < 300 pc, irrespective of the functional form we assume for K z . This is also among 
the first studies to find significant (> 2a) constraints on the local dark matter density, 
p DM = 0.0065 ± 0.0023 M pc" 3 (0.25 ± 0.09 GeVcm" 3 ). Fig. M shows the expected 
degeneracy between total surface density of visible matter and dark matter density: 
lower corresponds to a higher pdm for a given total surface density. However, these 
uncertainties do not include the systematic errors from misestimated photometric distances. 
For example, binary c ontamination w ould make the inferred distances smaller, and the 



distance estimation of 



An et al 



smaller than the one of 



Ivezic et al. 



(200 9j) wh ich we adopted in the present work is also ~ 9% 



(120081 ). We explore distance systematic in particular 
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for their implications on the dark matter density: if a 10% distance systematic error is 
considered, the infered change of /?dm is only 10%. 

Overall, our results are consistent with earlier studies, as Fig. [11] shows, which 
compares our inference about K z and about the total surface density E tot (z) with 
previous studies. The value derived by KG91 is S to t,|«|<i.i k pc = 71 ± 6 M pc~ 2 and 
p DM = 0.01 ± 0.005 M pc~ 3 . In the calculation of HF04, £ to t,| 2 |<i.i k pc = 74 ± 6 M pc~ 2 
and pdm = 0.007 M pc~ 3 , whereas the latter results are in a good agreement with ours. 

However, G 12 have re-analyzed the K dwarf sample from KG89b using the 



Jeans equation (IGarbari et al.l l201ll ) and they found a much large value of pdm 
0. 025^q;q}3 M q pc~ 2 (0.85to[% GeVcm -3 ). G12 use similar data as we do in this paper, 
but the sample is approximately five times smaller and the data have neither [Fe/H] nor 
[a/Fe] measured; the KG89b sample, however, has a simpler spatial selection function 
than the SEGUE K dwarf sample. There are a number of reasons why G12 may be 
overestimating the local dark-matter density: as they do not have measurements of [Fe/H] 
for the stars in their sample, their photometric distances are less accurate and less precise, 
and the metallicity distribution that they adopt from SEGUE has not been corrected for 
selection effects. The corrected metallicity distribution has a vertic al gradient close to the 



— 0.3dexkpc 1 gradient assumed by KG91 ( iSchlesinger et al 



20121 ): Appendix A of G 12 



shows that this leads to a slightly lower pdm- G12 also use a strong prior for the structure 
of the baryonic disk through their use of 15 isothermal baryonic components with fixed 
velocity dispersions. This model requires "dark matter" to account for any imperfections in 
the disk model. Most importantly, G12 cannot separate sub-populations based on elemental 
abundances as we do here, such that they may potentially mistake population gradients in 
a z (z) and h z (z) for the signature of dark matter. 

On the other hand, our estimates are in very good agreement (Fig. [T2|) with the recent 
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determination of pdm by BT12, who re-analyzed the data by iMoni Bidin et al.l ( 20 



They used the assumption of dV c /dR = to correct the model of 



Moni Bidin et al 



2L 



((2012) 



and derived a dark matter density of p DM = 0.008 ± 0.003 M pc~ 3 (0.30 ±0.11 GeVcm~ 3 ). 

Fig. [12] puts the recent local pdm estimate from K z constraints into context. The 
histogram with the gray shading in its center represents the results from G12, the one with 
black is from BT12 and the one with red is from the present work (hereafter Z12). The 
BT12 histogram reflects their pdm PDF including systematic errors; for Z12 several of these 
error sources have been addressed systematically, and the distance uncertainties have been 
incorporated as a systematic error term in the Z12 histogram in Fig. [T2J Therefore, the 
results shown in Fig. [12] from Z12 and BT12 should be on comparable footing. The joint 
PDF emerging from these three recent experiments is shown as the thick black histogram 
in Fig [121 indicating p DM = 0.0075 ± 0.0021 M p C - 3 (0.28 ± 0.08 GeVcm" 3 ). 

One can put these 'local' estimates of the dark matter density into the context of the 
expectations from the global fits of the Galactic dark matter halo: with a sp herical NFW 



cold d ark matter density profile (jNavarro et al. 



1996) and the parameters of 



Xueetal 



(120081) (virial mass M vir = 0.91±g^| x 10 12 M Q , virial radius r vir = 267l 2 j kpc, and c = 12.0), 
one would expect that the dark matter density at R Q is Pdm(-Rq) = 0.0063 M Q pc~ 3 , which 
is in very good, perhaps even fortuitous, agreement with our present work. 

Taken together with the existing work, our results continue to point towards a picture 
that the local Galactic potential is dominated by a fairly thin layer of stars and gas plus a 
dark matter density that is in accord with global fits to Galactic halo. 

It may seem surprising, why our limits we derive are not much tighter than those 
obtained by KG89b and subsequent work, based on smaller samples. This is in good part 
due to the fact that we used far fewer prior constraints on the models. E.g. we fit for 
the mass scale-height of the disk (zh) rather than assume it; similarly, we fit for the scale 
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height of the tracers simultaneously with the kinematic fit; and we have explored a range of 
functional forms for K z . Marginalizing over these factors, apparently leads to a seemingly 
similar error on the parameters than other more restricted analysis with smaller samples. 

6. Conclusion 

We have analyzed the K-dwarfs from SDSS/SEGUE for determining the total surface 
mass density and dark matter density in the local Galactic disc. At first, we divide our 
sample into three sub-populations through their [Fe/H] and [ct/Fe]. After considering 
the spatial selection function, the Galactic vertical number density profiles for different 
sub-populations are inferred from star counting. Secondly, we use maximum likelihood 
and Gauss-Hermite series to calculate the vertical velocity dispersion profile of each 
sub-population. Then different parameterized K z forms are used to solve the 'vertical' 
Jeans equation. We fit the observed vertical number density and vertical velocity dispersion 
profiles of the three sub-populations simultaneously, using MCMC to recover the PDFs 
of the parameters of the 'vertical force law'. In our results, for each sub-population, 
the vertical number density is approxmately single exponential profile and the vertical 
velocity dispersion is nearly isothermal. The scale height of number density profile and 
the vertical velocity dispersion increase with the decline of metallicity. Presuming that 
there is a thin gas layer with E gas = 13 M Q pc -2 , we derive a total surface mass density 
of 67 ± 6 M Q pc~ 2 at 1.0 kpc from the mid-plane, of which the contribution of all stars is 
42 ± 5M pc~ 2 and we infer a local dark matter density of pdm = 0.0065 ± 0.0023 M Q /pc~ 3 
(0.25 ±0.09 GeVcm- 3 ). 
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A. Other Model Functional Forms for K z 

Besides the two models for K z described in § 13.31 we also explore other three models 
to determine the total surface mass density and the mass density of dark matter, with the 
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purpose of checking whether the astrophysical inferences depend on the choice of these 
functional forms: 



hi) Error function model 



K z = 2itG 



erf 



z h 



ftas 



+ 4 71 G p DM z 



(Al) 



iv) General model I This form is an extension of KG model. We assume K z as the 
form of 



K z = 2ttG 



gas 



+ 4 7T G pbm z 



(A2) 



v) General Model II To reduce the number of free parameters, we set j3 = 7 and get 



K 7 = 2ttG 



+ S 



gas 



+ 4 n G pbm z 



(A3) 



We followed the same parameter estimation approach for these models and calculated 
the total x 2 between observed a 2i0 b s and v+ and the model predictions. The error function 
model and the general model II yield similar results for the surface density of visible matter 
and dark matter density as the KG89 and exponetial models discussed in the main text. 
Only the general models I, where the exponents 7 and j3 are independent yields to a degree 
of parameter degeneracy that makes the resulting K z {z) difficult to interpret. But overall 
this confirms that among physically plausible families of K z , the particular choice matters 
little. 
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Fig. 1. — The unweighted number distribution of the 9157 K-dwarfs from SDSS/SEGUE 



used in our anal ysis, shown in t 



from the SSPP ( Lee et al. 



re [a/Fe]- [Fe/H] abundance space, with abundances taken 



2008a 



bj). As volume completeness corrections have not yet been 



taken into account, this distribution overem phasizes the metal-poor, a-enhanced stars (see 



Liu fe van de Yen 



2012 



Bovy et al 



2012al ). Black boxes indicate the three abundance- 



selected sub-populations. I: metal-rich ([Fe/H] G [-0.5, 0.3], [a/Fe] G [0., 0.15]); II: in- 
termediate metallicity ([Fe/H] 6 [-1.0, -0.3], [a/Fe] G [0.15, 0.25]); III metal-poor ([Fe/H] G 
[-1.5, -0.5], [a/Fe] G [0.25, 0.50]). 
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bright plates faint plates 




Fig. 2. — Spectroscopic success rate for K-dwarfs in SDSS/SEGUE, used to estimate the 
selection function (§ I3.1.2p . The panels show the comparison of the distribution of the 
K-dwarfs that have good spectra from SEGUE (solid line) with the distribution of the 
presumably complete distribution of all photometrically detected point-source in the plate 
that satisfy the color-magnitude selection criteria. 
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Fig. 3. — Derivation of the spatial tracer number density for each abundance-selected sub- 
sample. From top to bottom, plots show the metal-rich, intermediate metallicity, and metal- 
poor sub-populations, respectively. From left to right, the panels show the actually detected 
number of stars as a function of galactocentric distance and height above the plane, the 
effective survey volume, and the implied number density (see § I3.ip . 
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Fig. 4. — Example fits of the vertical velocity dispersion, shown for three abundance selected 
samples in the bin 600 pc< \z\ < 700 pc. For all panels, the histograms show the observed 
vertical velocity distributions, the red solid curves are Gaussian fits without considering 
observed errors, and the green and blue solid curves are three and four moments Gauss- 
Hermite accounting for the observational errors. Without taking the individual errors into 
account overestimates the vertical velocity dispersion by about 2-4 km/s. 
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Fig. 5. — The derived vertical density profiles for the three abundance-selected sub-samples, 
with the model predictions from two of our parameterized models for the of K z force profiles 
(the KG89 and exponential model, respectively; see Eq. [ITl &; ITS]) . Filled circles are values 
estimated directly from the observations, dashed lines are model predicted values and shad- 
ows are the 68.3% errors in the recovered value of parameters. Red, green, blue symbols 
correspond metal-rich, intermediate metallicity, and metal-poor sub-population, respectively. 
xl is calculated from the observed z/* j0bs and modeled z/^ mod (Eq. |2^|) . 
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Fig. 6.— 



Analogous Figure to Fig. [51 but for vertical velocity dispersions. 



-37- 



KG89 Model 



Exponential Model 




300 600 900 
IZI [pc] 



1200 



300 600 900 
IZI [pc] 



1200 



Fig. 7. — The vertical dependence of the K z force of our two fiducial models, which include 
the stellar disk, a (thin) gaseous disk and a dark matter halo term. Left is the KG89 model 
(Eq. H7|) and right is the exponential model (Eq. IT8]) . In both panels, the dashed fat line 
shows K z (z) for the most likely parameters, and the band of grey points show a la sampling 
of the PDF for K z . 
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Fig. 8. — The PDF of model parameters and Zh for the models in Eq. [17] & [T8j On the 
left, the parameters for the KG89 model, on the right for the exponential model. Yellow, 
blue, and purple shades are 68%, 95%, and 99% confidence region, gray histograms are the 
marginalized PDFs for the individual parameters, and red lines represent the most likely 
value of each parameters. 




Fig. 9. — Analogous Figure to Fig. El but for the PDFs of and the fitted scale height of 
the tracer population of stars, h z , illustrated for the case of the metal-rich sub-population. 
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Fig. 10- 



Analogous Figure to Fig. [HJ but for the PDFs of pdm and E*. 
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Fig. 11. — Comparison of the vertical force K z (z) and corresponding surface mass density 
2tot.<|z| implied by the best fits of the various model families to our data. The dot ted, the 
dark solid, and the long dashed red lines represent three cases of our predictions based on 
KG89 model family: the one without DM (as a limiting case), the one with best fitting DM, 
and the one with pdm = O.OO8M pc -3 (from BT12), respectively. The 68% uncertainty 
intervals on the surface- mass density are shown at a few representative points. The \z\ range 
where the lines are drawn thicker represent the \z\ range of our sample stars. Dashed line is 
the same model but from the prediction of KG89. 
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Fig. 12. — Comparison of normalized PDFs of pdm derived by three independent works. The 
dashed dot line and gray shade represent the results of G12, dashed curve and shade are 
from the calculation of BT12, and red curve and shade are the results of our present work. 
Thicker black line means the joint PDF emerging from these three histogram. 
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Table 1: Fitting results of the two models, with E gas = 13 M Q pc 2 (IFlynn et al.ll2006l ) 



KG89 model (Eq. d7j) Exponential model (Eq. M 

S ,star [M pci 42 ±6 41 ±5 

Pdm [M pc- 3 ] 0.0064 ± 0.0023 0.0060 ± 0.0020 

z h [pc] 2451*1 200^ 

h zA a [pc] 259 ± 12 260 ± 15 

h z , 2 a [pc] 450 ± 26 465 ± 33 

h z , z a [pc] 852 ± 30 910 ± 71 

cr ,i a [km/s] 15.4 ± 1.3 15.8 ± 1.3 

a ,2 a [km/s] 23.0 ± 2.0 23.6 ± 2.0 

a ,3 a [km/s] 34.2 ± 2.2 35.8 ± 2.2 

xl 25.23 26.36 



X 



I 15.80 14.54 



L 41.03 40.90 



Xtot 



"Number 1, 2, and 3 represent metal-rich, intermediate metallicity, and metal-poor sub-sample, respectively. 



